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Abstract 

We present a new determination of the QCD strong coupling constant 
based on precise lattice calculations of the T spectrum. The largest sys- 
tematic uncertainty in previous such determinations resulted from the 
absence of vacuum polarization from light quarks. We substantially re- 
duce this error by including two flavors of dynamical light quarks and 
extrapolating to three. We find a^'(8.2 GeV) = 0.196(3) for three light 
flavors, corresponding to a^(Mz) — 0.115(2). This is significantly more 
accurate than previous determinations using this or any other technique. 

PACS numbers: 12.38. Gc, 12. 38. Aw, 12.38.Bx, 14.40. Gx 

In recent years it has become apparent that numerical simulations of lat- 
tice quantum chromodynamics (QCD) could provide accurate and reliable de- 
terminations of the strong interaction coupling constant |l| ^ ^, |^, |j. The 
largest systematic error in previous lattice determinations resulted from the ne- 
glect of quark vacuum polarization, whose sizable contribution was estimated 
perturbatively. In this paper we significantly reduce this systematic error by 
presenting simulation results that include dynamical light quarks. We find 
av'(8.2 GeV) = 0.196(3) for the strong coupling constant with = 3 flavors 
of light quarks; V denotes the physical scheme based on the static-quark po- 
tential, discussed in Refs. [|[ 0]. The corresponding MS coupling at the Z mass 
is a^{Mz) = 0.115(2), consistent with the world average of 0.117(5) §, but 
significantly more accurate. 
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In this study we used two different sets of simulations, one without quark vac- 
uum polarization (n/ — 0), the other with two flavors of light quarks {nj — 2) 
The nj — 2 gauge-field configurations were generated using a staggered-quark 
action for light quarks with the Hybrid Molecular Dynamics algorithm, and 
both used the standard Wilson action for gluons. For b quarks, we employed 
the nonrelativistic formulation of quark dynamics (NRQCD); details are in 
Refs. 

There are two steps in a determination of the coupling constant from lattice 
simulations. The first is to specify or accurately determine the parameters of 
the lattice lagrangian, and in particular the lattice spacing a. The second is to 
use this lagrangian to compute nonperturbatively an appropriate short-distance 
quantity. Comparison with the perturbative expansion for the same quantity 
fixes the coupling. 

Just as for continuum QCD, the bare coupling constant and masses must be 
provided as input. In the lattice action, the bare coupling yiat appears in the 
parameter /3 — The lattice spacing is not an input. Rather, it is specified 

implicitly by (3; for each (3 there is a corresponding a. By expressing all dimen- 
sionful quantities in units of a, it is scaled out of the action, and so serves to set 
the overall mass scale. As a result, the simulation produces a value for some par- 
ticular mass M only in the dimensionless combination aAI. We must therefore 
know a before we can compare M to its experimental value. In our simulation, 
we compute the mass difference between the T and T' mesons, aAAf(T' — T), 
and between the T and the spin average of the Xb states, aAM{xb ~ T). We 
then divide these by the experimentally measured mass differences to obtain 
two independent estimates for a. 

Heavy-quark systems possess several properties which permit us to measure 
a accurately Q. They are essentially nonrelativistic; the use of a nonrelativistic 
effective action to exploit this allows a large portion of the spectrum to be 
computed efficiently and precisely. They are physically small, and do not suffer 
from finite-volume errors on modestly sized lattices. Their spin-averaged mass 
splittings are observed experimentally to be nearly independent of the heavy- 
quark mass, varying by only a few per cent between the T and -0, making our 
results insensitive to tuning errors in the bare 6-quark mass. Because including 
vacuum polarization in simulations from nearly massless quarks is difficult, it is 
common to use unrealistically large u- and d-quark masses, and then extrapolate 
to the correct values. However, due to the small size of the T and the large 
momentum transfers between its constituents, the bare light-quark masses of 
about 25 MeV in our Uf — 2 simulations are negligible, making extrapolation 
unnecessary. Finally, T decay rates are negligible as compared to their energy 
splittings, making the effect of light-quark mass values on nonanalytic threshold 
behavior unimportant. 

Our results for nf — and 2 are summarized in Table ^ The 6-quarks 
in T, T' and Xb mesons typically exchange momenta of order 1 GcV, so that 
the appropriate number of light flavors to include in a study of their dynamics 
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1.80 


.237(10) 


2.38(10) 


.178(5) 


2.47(7) 



Table 1: Lattice QCD simulation results for the difference between the T' and 
T masses, and between the spin-averaged Xb mass and the T mass. Results are 
given for different bare gluon couplings /3 — Q/gf^^ and bare quark masses M°, 
and for Uf = and 2 flavors of light quarks. The correct bare mass for the 
6-quark is aAI^ — 1.7(1). The splittings are corrected for 0{a^) errors in the 
gluon action. The errors shown are statistical and result from our use of Monte 
Carlo methods in the simulations. Values for a^^ are in GeV and are obtained 
using AM(T' - T) = 0.563 GeV and AM{xb - T) = 0.440 GeV. 

is three. Having results for n/ of both and 2 will allow us to accurately 
extrapolate to Uf ~ 3. In addition, we varied the bare &-quark mass around 
the correct value of aM° — 1.7(1). As expected, the splittings showed little 
sensitivity to Af°. 

We have determined both by fitting the splittings separately and simul- 
taneously. Fitting separately produced a discrepancy in of a couple standard 
deviations for nf = 0, and about half this for nf = 2. As we will show, this small 
discrepancy vanishes when we extrapolate n/ to three. Insofar as the effect is 
real, it is likely due to the larger intrinsic momentum transfers for S states as 
opposed to P states.[] Fitting simultaneously, we find that our simulation data 
are consistent with w 2.4 GeV. Our results for other low- lying excitations 
and spin splittings of the T system, using this value for a~^, are displayed in 
Figures |l| and ^ where they are compared with their experimental values. The 
excellent agreement supports the reliability of our simulations. We emphasize 
that these are calculations from first principles; our approximations can be sys- 
tematically improved. The only inputs are the lagrangians describing gluons 
and quarks, and the only parameters are the bare coupling constant and quark 

^ Perturbation theory, though not justified at these momenta, provides a quafitative expla- 
nation of the effect on the determination of a~^. The centrifugal barrier makes the average 
separation between the quarks in the P-state Xb larger than for the S-state T or T', as is fa- 
miliar from hydrogen or positronium. Consequently, the typical exchanged momentum for Xb 
quarks, q^i, , is smaller than q-f' ■ The perturbative binding energy is given by (q)C'p Mi, /16, 
with q = qyi for T' and q^j^ for Xb- Since q^^ < Qf 'i the Xb is more tightly bound. However, 
for = 0, this effect is exaggerated, as o^'(g) increases more quickly than Q^'(g) with de- 
creasing q. Thus, for < 3, AM(x6 — T) should be underestimated relative to AM(T' — T), 
as is observed. Fitting to data would then require a larger for AM(x(, — T) than for 
AM(T' - T). 
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Figure 1: NRQCD simulation results for the spectrum of the T system, including 
radial excitations. Dashed lines indicate experimental values for the triplet S- 
states, and for the spin-average of the triplet P-states. The energy zero from 
simulation results is adjusted to give the correct mass to the T(l'^S'i). Results 
are from a simulation with Uf — Q (filled circles) and from one with Uf — 2 
(open circles), using — 2.4 GeV for both. The errors shown are statistical; 
systematic errors are several tens of MeV. 



mass. In particular, these simulations are not based on a phenomenological 
quark potential model. 

The 5-quark action includes relativistic corrections in the quark velocity v to 
0{v'^), and corrections to ©(a^) in the lattice spacing. The leading systematic 
error affecting is almost certainly due to finite-lattice-spacing errors in the 
gluon action. To correct for this requires an additional 0{a^) interaction. As 
it is sensitive to short distances, the effect of this interaction can be estimated 
using perturbation theory, which gives a mass shift of 

Here qs ~ 5 GeV is the typical momentum transferred through the interaction, 
and V'(O) is the meson wavefunction evaluated at the origin. Using Uf — 
simulation results for ip{r) and we find that aAMg is 0.0036 for the T 
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Figure 2: NRQCD simulation results for the spin structure of the lowest-lying 
P-states. Dashed lines indicate experimental values for the triplet P-states; the 
dotted line shows their spin average. Masses are relative to the spin-averaged 
state. Results are from a simulation with Uf — (filled circles) and from one 
with Uf — 2 (open circles), using a^^ = 2.4 GeV for both. The errors shown 
are statistical; systematic errors arc within about 5 MeV. 



and 0.0023 for the T' when aM^ — 1.71; there is no shift for the P-state x's. 
The correction increases approximately linearly with mass and is about 30% 
larger for nf — 2. The splittings in Table H include this correction (although 
Figures |l| and ^ do not). It shifts by almost a standard deviation in the 
case of the Xt — "^^ splittings, but is negligible for the T' — T splittings. 

To check the validity of Eq. (1), we examined the perturbative prediction 
for the hyperfine splitting between the T and ?7f,, given by the similar formula 



AAfM. = ^^^|^(0)r, (2) 



where Mi, is the pole mass of the 6-quark |Tl| . This formula gives a AMhfs = 
0.0122(10) when aM" = 1.71 and nj = 0, which compares very well with 
the nonperturbative result 0.0123(2) we obtain from our simulation jl^. This 
suggests that our perturbative estimates of the 0{a^) shifts (Eq. (1)) are quite 
reliable. We have also computed these shifts nonperturbatively using a lattice 
potential model, and found essentially identical corrections. Other systematic 
errors are higher order in a or in the quark velocity v, and so are most likely 
negligible. The dominant errors in determining a^^ are statistical. 

Having determined the lattice spacing, the second step is to use simulation 
results to extract the coupling constant. Our approach closely parallels deter- 
minations based on high-energy phenomenology, with the simulation playing 
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3.41 
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0.5708(1) 
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1.5123(2) 
1.8337(3) 
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.181 
.181 
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.1785 



Table 2: Perturbative and simulation results for several small Wilson loops [15]. 
The entries in the final two columns list the values for ay extracted by comparing 
the simulation results for each — In to its perturbative expansion, correct to 
second and third order, respectively. 

the role of an experimental measurement. We identify short-distance quantities 
whose perturbative expansions are known at least through second order, and 
then determine these quantities nonpcrturbatively using our simulation. By 
equating the perturbative expansion to the nonperturbative value, we can solve 
for the strong coupling constant. 

One of the most ultraviolet quantities in lattice QCD is the expectation 
value of the 1x1 Wilson loop operator. Its perturbative expansion to 0{a^) 
is [H 

-InM^i,! = y al"^^(3.41/a) {l - (1.185 + 0.070 n/) al"^^}. (3) 

The strong coupling constant was defined in Refs. ||, in terms of the 
static-quark potential. To simplify the discussion of larger Wilson loops that 
follows, it is convenient to regard Eq. (^), with no higher-order terms in ay, as 
defining av The two definitions differ only at C'(a^). The scale 3.41/a follows 
from the technique described in Ref . , and indicates the important momentum 
scale in Wi^i. It corresponds to 8 - 9 GeV for our lattices, confirming that Wi^i 
is very ultraviolet. 

Simulation results for — InW^i i are listed in Table |^. From Eq. (^) and the 
extracted from the Xb ^ ^ splitting we find 

a<,"> (8.76(41) GeV) = 0.1517, (4) 
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(n 



aV" (8.42(24) GeV) = 0.1785 (5) 

or equivalently, 

a, (8.2 GeV) - | o.l800(16) for = 2 , 

where the errors are due to the statistical errors in a~^. As mentioned above, 
we must extrapolate our results to n/ = 3. Perturbation theory suggests that 
l/av"^"* is more nearly linear for small changes in Uf than a^v^^\ and so we 
extrapolate the inverse couplings to obtain 

a^^^ (8.2 GeV) = 0.1959(34) . {xb - T} (7) 

Since we are extrapolating by only 9%, extrapolation errors are probably smaller 
than the statistical errors quoted here; perturbation theory indicates that they 
are about 0.2%, which is neghgible. 

Repeating this analysis using from T' — T gives 

'^Vs^GeV^- /^•1^04(22) forn^-0 , . 

(8.2 GeV) - I 0^779(23) for nj = 2 , 

and 

a^^' (8.2 GeV) = 0.1958(46) , {T' - T} (9) 

which agrees with the Xb — '^ value. We have then, from both Xb — "^ and T' — T, 
our primary result :0 

(8.2 GeV) = 0.1959(34) . (10) 

Because the internal momenta transferred by the T constituents are small 
relative to c- and 6-quark masses, it would be incorrect to extrapolate ay ob- 
tained with n/ = and 2 light quarks directly to = 4 or 5. The correct way 
to incorporate these heavier flavors is to run the coupling to below the c-quark 
threshold, then apply matching conditions as a is run back up through the c 
and b thresholds. We will use this procedure to obtain at the Z mass. 

We chose the 1x1 Wilson loop because we expected nonperturbative effects 
to be very small due to the large momentum scale it probes. To verify this, we 
have determined a^"'''' using larger Wilson loops and loops with different shapes, 
which should have significantly larger nonperturbative contributions. In Table ^ 
we give the perturbative expansion coefhcients for the smallest Wilson loops. 
These coefficients are defined by 

- In W^%^ - ^1"'^ [^v'\qn,m)] ' , (11) 



^We expect these two determinations to be statistically correlated, and quote as an error the 
uncertainty in the Xb ^ determination, rather than combining the errors as if independent. 
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0.151(7) 
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0.4884 
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0.1381 
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0.5214 


2.57(12) 


8.8(4) 


0.1517 
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Table 3: Results for a^"'(8.2 GeV) with jS = 6.4, 6.2 and 6.0 using the Xb - T 
splitting. The scale qn = 3.41/a; qn and a^^ are in GeV. 

with at"^^ defined by Eq. (^) We also quote simulation results for these 

quantities, and values for ai"^'' (3.41 / a) obtained by matching second- and third- 
order perturbative expansions to the simulation results. The n/ = results, 
when third-order perturbation theory is used, show that the small loops all 
give the same value for a^' (3.41/a) to within less than one per cent. This 
confirms that nonperturbative effects are completely negligible. (In fact, the 
slight variation is likely due to fourth-order perturbative corrections.) The rif = 
2 results are consistent with this conclusion, although the test is somewhat 
less stringent since the third-order perturbative coefhcients are not known for 
Uf ^ 0. 

To verify that additional lattice-spacing errors are under control, in Table ^ 
we compare our result for a'"' (8.2 GeV) computed from a simulation at /? = 6.0 
with those from /3 = 6.2 and 6.4 [|l^. The scales for the corresponding 1x1 
Wilson loops range from 8 to 14 GeV. That these give consistent values for 
a'"' (8.2 GeV) indicates that, within errors, the coupling constant is scaling 
correctly, and confirms results found in Ref. ||^. 

Eq. ( p^ is our final result. However, to facilitate comparison with other 
determinations we convert our result to the MS scheme, which is related to the 
V scheme by 



a^\Q)^ai^'\e'^'Q) {l + 2 a^'/^ + • (12) 



MS 



Our result is then equivalent to 

a|l(3.56 GeV) = 0.2203(84) , (13) 

with the error now dominated by the unknown third-order contribution to 
Eq. (|12|), which we estimate as (a^?')^ = 0.0075. 

We numerically integrated the third-order perturbative beta function for 
a^'^ and applied appropriate matching conditions at quark thresholds to evolve 
it to several other scales M: 



a 



(»/) 



0.304(17) for Q = 1.7 GeV « Mc and n/ = 3,4 



MS 



f>{Q) ^ { 0.203{7) for g = 5.0 GeV w Mb and n/ = 4,5 (14) 



0.115(2) for Q = 91.2 GeV = Mz and Uf = 5 
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Source 



Uncertainty 



Converting from a'y'' to agl 1.7% 

Statistical error in determination of .9% 

Extrapolation in n/ .2% 

Finite a and 0{v^) errors .2% 

Fourth-order evolution of aj7g .01% 



Table 4: Sources of error in a^{Mz)- 



The last of these results is consistent with, and significantly more accurate than, 
the current world average a^{Mz) = 0.117(5) . 

We believe we have accurately estimated sources of error. We have checked 
our result using two different mass splittings to determine a~^, and four inde- 
pendent Wilson loops to extract ai/(3.4/a). We have checked for consistency 
against variations in the quark mass and lattice spacing. Our estimate of the 
third-order perturbative contribution to the relation between and a-^ is 
consistent with first- and second-order contributions, and with the third-order 
terms in Table |[ Our result also agrees well with earlier determinations based 
on lattice simulations [Q, ^. These were performed without dynamical light 
quarks. To estimate the effect of light quarks, av\ computed at 3.41/a, was 
run down to a scale typical of momenta exchanged in the T or ^, determined as 
in Ref. |0, There it was equated to a^', since it is at this scale that the two 
are required to produce the same values for splittings, and av' was then run 
back up to the desired scale. We may repeat this procedure as an alternative 
to a direct extrapolation in ny, now using results for both zero and two light 
flavors. This gives a^(M2) = .112(4) from n/ = data, and .115(3) from 
Uf — 2. So this method yields values consistent with, but not as reliable as, 
extrapolation in ti/ . 

There are prospects for substantially improving the accuracy of our result 
fairly soon. Sources of error in our value for a^{Mz) are listed in Table ^. The 
dominant error is unrelated to our primary result of Eq. (10), but rather is due 
to the conversion to MS. The total error could be cut in half by computing the 
third-order correction to Eq. (^2|), a straightforward perturbative calculation. 
The error in will decrease with improved statistics; we have now completed 
a new simulation that should soon reduce it to about .6%. Use of an improved 
gluon action would remove the need for the correction in the xti — T analysis, 
at little additional cost jl^. Finally, a simulation with either nj — 3 or A light 
quarks would eliminate the extrapolation error and would require roughly the 
same amount of time as for n/ = 2. 

There are a variety of additional lattice calculations that would provide a 
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broader check on the consistency of our result. The first would be to repeat 
our determination using charmonium, extending previous studies by including 
dynamical light quarks. A similar analysis is also possible using light hadrons. 
For reasons outlined in the introduction, systematic errors for light hadrons are 
not nearly as well understood as for heavy mesons. It is important that control 
of these errors be improved to demonstrate that these systems yield consistent 
results. Finally, simulations which include light-quark vacuum polarization are 
still relatively rare. It would be very useful to repeat this analysis using an n^^ 
other than two, algorithms for generating gauge-field configurations other than 
HMD, and Wilson fermions rather than staggered. 

In this paper we have demonstrated that lattice simulations provide one of 
the simplest, most accurate, and most reliable determinations of the strong cou- 
pling constant. The fact that a lattice simulation of nonperturbative hadronic 
structure at scales smaller than 1 GeV agrees with perturbative analyses of 
high-energy jet formation is striking confirmation that these diverse phenomena 
are governed by a single theory — QCD. Furthermore, this result demonstrates 
our growing mastery over both the nonperturbative and perturbative aspects of 
the theory. 

An independent lattice determination of a^{Mz) with dynamical light 
quarks has recently appeared in Ref. j2^. They obtain results consistent with 
ours, though with larger errors. 

We thank Aida El-Khadra, Paul Mackenzie, Urs Heller and Jochen Fingberg 
for several useful discussions; Greg Kilcup and collaborators, and the HEMCGC 
collaboration, for generously providing us with their gauge-field configurations; 
and Jonas Berlin for help adapting the HEMCGC configurations for the Cray. 
Our simulations were carried out at the Ohio Supercomputer Center and the 
National Energy Research Supercomputer Center. This work was supported by 
grants from the NSF, the SERC, and the DOE (DE-FG02-91ER40690, -FC05- 
85ER250000, -FG05-92ER40742, -FG05-92ER40722). 
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